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Abstract 


The deterministic particle transport code HZETRN was developed to enable fast and 
accurate space radiation transport through materials. As more complex transport 
solutions are implemented for neutrons, light ions (Z < 2), mesons, and leptons, it is 
important to maintain overall computational efficiency. In this work, the heavy ion (Z > 

2) transport algorithm in HZETRN is reviewed, and a simple modification is shown to 
provide an approximate 5x decrease in execution time for galactic cosmic ray transport. 
Convergence tests and other comparisons are carried out to verify that numerical 
accuracy is maintained in the new algorithm. 

1. Introduction 

The deterministic particle transport code, HZETRN (High charge (Z) and Energy TRaNsport) 
[Wilson et al. 1991; Slaba et al. 2010a, 2010b], is a fast and accurate computational tool for transporting 
space radiation through shielding materials and human tissue. Recent updates to the code have provided 
significant reductions in run time [Slaba et al. 2010a], and further work reduced overall discretization error 
to less than 4% for shielding materials and thicknesses commonly found in space applications [Slaba et al. 
2012a]. This increased efficiency has enabled new calculations that would have otherwise been 
computationally prohibitive. First, a ray-by-ray transport methodology has been developed wherein 
transport is explicitly performed along each ray direction in a shielding thickness distribution, and no 
constraint is placed on the number or ordering of materials along a ray [Slaba et al. 2011a]. Second, the 
transport of pions, muons, electrons, positrons, and photons has been added into HZETRN at a moderate 
computational cost [Norman et al. 2012a, 2012b; Slaba et al. 2012b]. Third, calculations may be performed 
through detailed vehicle geometry on a minute-by-minute basis along a trajectory. This type of analysis has 
been useful in providing more rigorous validation comparisons and uncertainty estimates [Slaba et al. 
2011b, 2012b], 

Further developments that lead to more efficient numerical methods or algorithms are beneficial 
and may also enable new analyses to be completed. In this work, the heavy ion (Z > 2) transport algorithm 
in HZETRN is reviewed, and a simple modification is made that allows interpolation to be avoided in the 
ion collision source term. It is shown that the modification provides an approximate 5x decrease in 
execution time without sacrificing numerical accuracy. 

2. Heavy Ion Transport Formalism 

In this section, the heavy ion transport formalism in HZETRN is reviewed. The Boltzmann 
transport equation with the continuous slowing down and straight ahead approximations is given as 
[Wilson etal. 1991] 


B^E)} = x; J7 
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In equations (1) - (3), (f>.(x,E) is the flux of type j particles at depth x with kinetic energy E, Aj is the 

atomic mass number of a type j particle, Sj(E) is the stopping power of a type j ion with kinetic energy E, 
cij(E) is the total macroscopic cross section for a type j particle with kinetic energy E, and Ujk(E,E') is the 
differential macroscopic production cross section for interactions in which a type k particle with kinetic 
energy E produces a type j particle with kinetic energy E. The summation limits in equation (1) are taken 
over all projectiles that might produce a given fragment. The boundary condition spectrum, / .(E) , is 

considered to be a known function over a broad energy spectrum. 

The ion stopping powers, S/E), in equation (2) may be approximated in terms of the proton 
stopping power, S{E), as [Bethe et al. 1930] 

j-SjWnvjSiE), (4) 

where the scaling parameter is expressed in terms of the ion mass and charge (Aj, Z ; ) according to 
v j = Z 2 j / Aj . The transport equation may now be written as 

B[^ (x, E)]=J2 J e Vjk (£> E 'Mfc (*> E ') dE ' ' , (5) 

k 

with the modified linear differential operator 

B[^{x,E)] = A S (E ) + a j (E) . (6) 

Approximating the ion stopping powers with the scaled proton stopping power allows simplified numerical 
algorithms, as has been discussed elsewhere [Wilson et al. 2006]. 

For heavy ions, it is noted that projectile fragments have an energy and direction very near that of 
the projectile, while target fragments are produced nearly isotropically with low energy and travel only a 
short distance before being absorbed [Wilson et al. 2006]. The approximate decoupling of target and 
projectile fragments is discussed in detail by Wilson et al. [1995] and suggests that target fragments can be 
neglected in the heavy ion transport procedure (their contribution to dose is approximately accounted for 
after the transport procedure). The equal velocity assumption for heavy ions can be expressed in the 
production cross section as [Shinn et al. 1992] 

a jk (E,E') = a jk (E')8(E-E'), (7) 

where (/j k (E) is the partial fragmentation cross section for interactions in which a type k particle with kinetic 
energy E produces a type j particle. In this case, equation (5) becomes 

B[cf> j (x,E)\=^a jk (E) ( f> k (x,E). (8) 

k 

The absence of target fragments in the heavy ion transport procedure allows one to take the summation in 
equation (5) over all projectiles with mass greater than that of the fragment. If all the transported particles 
are ordered according to mass, then equation (5) can be succinctly written as 

B[^ (x, E)] = ^2 a jk (E) <j) k (x, E ) , (9) 

k>j 
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which is the transport equation found in Wilson et al. [2006] and will be referred to as the heavy ion 
transport equation. The upper summation limit in equation (9) can vary, but it is common to use no fewer 
than 59 ions [Cucinotta et al. 2006]. 

For light particles (n, p, 2 H, 2 H, 3 He, 4 He), both projectile and target fragments are included in the 
transport procedure. The broad energy distribution in collision events also indicates that the equal velocity 
assumption in equation (7) cannot be used. In this case, no simplifications to equation (5) are available, and 
the summation is taken over all light particles. It should be noted that for solar particle events with a 
negligible heavy ion component, only the solution to the light particle transport equation is required. For 
galactic cosmic ray (GCR) environments, there are both heavy ion and light ion components, and solutions 
to the light particle and heavy ion transport equations must be evaluated simultaneously. The light particle 
transport solutions and coupling to heavy ion transport has been discussed in detail elsewhere [Slaba et al. 
2010a] and will not be repeated here. 

3. Heavy Ion Marching Equation 


In this section, the mathematical formalism that allows equation (9) to be inverted and expressed 
as an explicit marching procedure is reviewed. The following derivation was taken from Wilson et al. 
[1977, 1986, 1991, 1995, 2006] and Shinn et al. [1992]. The scaled heavy ion flux is defined as 

j (z, r) = v j S{E)(j) j (x,E), (10) 

where r is the proton range defined in terms of the stopping power according to the equation 


_ r E dE' 

? = Jo S(E ') ' 

Equation (9) may be written in terms of the scaled flux and proton range as 


( 11 ) 
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where the quantities, &j(E) and cr- k {E) have been replaced with a } (r) and cr jk (r) since, for a given r, 

equation (11) may be inverted to find E. Equation (12) may be inverted using the method of characteristics 
to obtain the Volterra integral equation 

ipj(x,r) = + i>jx) 

+J2~f e ~^^ r ' X ^ a jk ( r + V j X - x' ,r + v j x ')dx' , 

k>j v k Jo 


with the exponential term 


0j(r,x) = f o cr.(r + vfidt . (14) 

The complete details of how equation (9) is transformed into equation (12) as well as the inversion of 
equation (12) have been given by Slaba et al. [2010c]. Equation (13) may be written as a marching 
procedure in terms of a step-size, h, as 
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(15) 


V’jO r + h,r) = e W r ’ h ' > i(> j (x,r + vfr) 

+J2~f jk (r + VjX')il> k (x + h - x\r + v^dx' . 

k>j u k J ° 

The integrand in equation (15) is approximated by considering the solution to the homogeneous form of 
equation (12). The homogeneous equation neglects secondary particle production through nuclear 
interactions and accounts only for the slowing down of particles due to atomic interactions and the loss of 
particles due to nuclear absorption. If the step-size is taken to be sufficiently small such that 


h « 




(16) 


(i.e. much less than the nuclear mean free path), then the local truncation error will be negligible as the 
particles will not have travelled far enough to participate in a nuclear collision [Wilson et al. 1991]. The 
homogeneous solution is given by the equation [Wilson et al. 2006] 


^j(x + h,r) = e l3j< ' r ' h) 'ip j (x,r + v^h), 

which may be equivalently written as 

ip k (x + h - x\r + v-x') = e A( r + V j x ,h x ^ k (x,r + VjX'+ v k (h - x’)) , 

and substituted into equation (15) to obtain 

' ipAx + h,r) — e~^^ r ’ h ' > tl)Ax,r + v-h) 


(17) 
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The fragmentation cross section and particle flux appearing in the integrand of equation (19) are 
approximated as 


(19) 


a jk( r + v j x ') « a A r +v i h / 2 )’ 


( 20 ) 


i’k( x ’ r + v j x '+ v ki h - x '))~ M x ’ r + + u k) h / 2 ) • 


( 21 ) 


Equation (20) is accurate since heavy ion fragmentation cross sections vary slowly as a function of energy 
(range), and most collisions occur at energies for which the ion range is large compared to typical step- 
sizes. Equation (21) is accurate since the scaled heavy ion fluxes are nearly constant at low range values. At 
higher range values, the step-size is comparatively small, and the step-size perturbation has a small affect. 
These statements are supported in Figure 1, showing representative fragmentation cross sections and scaled 
ion fluxes as a function of proton range (bottom axis) and kinetic energy (top axis). 
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Figure 1. The left pane shows the partial fragmentation cross section for 56 Fe + A1 as a function of proton 
range and kinetic energy. The right pane shows the scaled ion fluxes at 20 g/cm 2 in aluminum exposed to 
solar minimum GCR conditions given by the Badhwar-O'Neill 2010 model [O'Neill 2010] with a solar 


modulation parameter of 475 MV. 

Upon substitution of equations (20) and (21), the heavy ion transport equation becomes 


ipj(x + h,r) 


— e + v-h) 

+E— CT i*( r + v i h / 2 )M x > r + 

k> j V k 



e -0 j (r,x')-0 k (r+v j x\h-x') 


dx '. 


( 22 ) 


The remaining integral is evaluated by expanding the exponential terms in the integrand of equation (22). 
The complete details of the expansion and integral evaluation have been given by Slaba et al. [2010c]. The 
integral is approximately evaluated as 
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Jo 


-PJr^'yPkir+va'A-x 1 ) 


dx ' 


e ~a j {r+v j h/2)h _ e ~a k (r+(v j +i' k )h/2)h 
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The final marching procedure is now given by the equation 

ipj (x + h,r) — e l3 ^ r ’ h) 'i/} j (x,r + v-h) 

__ v • . x (24) 

+Y,— a jk ( r + v i h / 2 )A( x ’ r + E + v k) h / 2 ) c ' J x( r )> 

k> j v k 

This is the heavy ion marching equation given by Wilson et al. [2006]. The heavy ion fluxes are evaluated 
on a pre -defined energy grid that is converted to proton range values. There are typically 100 grid points 
spaced logarithmically from 0.01 MeV/n to 50 GeV/n. In order to obtain the ion flux at a grid point, r b 
interpolation is needed to evaluate the terms ip .(x,r + v .h) and ip k (x,r + {v . + v k )h / 2) . The 

interpolation method that has been typically used in FIZETRN is log-log cubic Lagrange [Wilson et al. 
1995]. In this method, the logarithm of both the dependent and independent variables is evaluated. A cubic 
Lagrange interpolation scheme is applied, and the resulting value is exponentiated. 
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4. New Numerical Method 


In this section, the original numerical implementation of the heavy ion marching equation is 
modified to provide a more efficient algorithm. As stated in the previous section, interpolation is required 
to evaluate the ion fluxes at values not on the prescribed energy/range grid. By examining the summation 
term in equation (24), it can be determined that the number of interpolations per marching step scales 
linearly with the number of energy grid points and quadratically with the number of isotopes. For typical 
calculations, this number is on the order of 10 3 interpolations per step. For a given range value, r,-, an 
alternative strategy to evaluate equation (24) is to use the nearest-neighbor approximation 

i>k( x ’ r i + + v k) h / 2 ) ~ i>k( x ’ r i) ’ ( 25 ) 

where r* is the nearest grid point value to r + (i/. + v k )h / 2. For each r u the heavy ion marching 
equation becomes 


Ip Ax + h,rA = e ^’^ipAx,^ + v h) 


+Y,— a jk{ r i+ u j h / 2 ) A (*> r i ) C jk ( r i )• 


(26) 


k>j v k 


In this equation, no interpolation is needed to evaluate the ion flux within the production 
summation. Consequently, the only interpolation required is for i/f (x. r + v.h ) , and the total number of 

interpolations per step is now directly proportional the number of energy grid points. In order to take 
advantage of the nearest neighbor approach, the mapping from each r t to the nearest neighbor of 
r t + (v j + v k )h / 2 is evaluated once and used repeatedly throughout code execution. As in the previous 

section, equation (25) is accurate at low range values since the scaled ion fluxes are nearly constant in this 
region. At higher range values, the step-size perturbation has a small effect. Further, the number of grid 
points typically used for heavy ion transport provides sufficient fidelity so that the nearest-neighbor range 
value is not far from the actual perturbed value. 

It should be noted that at the first step from the front boundary, secondary ion fluxes are changing 
rapidly, and the previous heavy ion solution is more accurate than equation (26). Similar behavior is 
observed near material interfaces, but the impact is less noticeable than at the front boundary. In the 
updated heavy ion algorithm, a hybrid approach has been adopted. For the first step into any material, 
including material interfaces, the original heavy ion marching equation (24) is evaluated. All other steps 
use the simplified form in equation (26). Log-log cubic Lagrange interpolation is still used to evaluate the 
term ip j (x, r l + z 'jh) . 


5. Flux Comparison 


In this section, a flux-based comparison between the original and new heavy ion transport 
algorithms is given. In this comparison, both algorithms were used to transport a solar minimum GCR 
environment through a 100 g/cm 2 aluminum slab and a 100 g/cm 2 water slab, and fluxes were saved at 
various depths in both slabs. The solar minimum GCR environment was calculated with the Badhwar- 
O'Neill 2010 GCR model [O'Neill 2010] using a solar modulation parameter of 475 MV. The relative 
difference between the results computed with both algorithms is computed as 


R A E v x k) = 


^\ E ^ k )~^f\E v x k ) 


(27) 
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where ^^(E^Xf.) and 6]^ (E i ,x k ) are the fluxes for each ion, j, at each energy, E h and depth, x k , 

computed with the original ( O ) and new (AO algorithms, respectively. It was found that the absolute value 
of the relative difference between the old and new algorithms never exceeded 5% over all energies, ions, 
and depths. In order to summarize these results, the following quantities were computed as a function of 
energy, 


R min( E i) = mi nR(E v x k ) 

J,k 


R nmx( E i) = max Rj( E i,x k ) 

J,k 

, N, N, 


R (E.) = 

avg 1 N-N 


-E£W^) 


(28) 

(29) 

(30) 


j x k= 1 j = 1 




N, JV,- 


EE[ R j( E v x k)- R , vg ( E u 


N j N x k^l^l 3 


(31) 


Equations (28) - (31) have also been computed as a function of particle type and depth. The results are 
shown in Figures 2-4. The symbol in each figure is the average of the relative differences, R mg (E,). The 
outer-most error bars represent the minimum and maximum relative differences defined in equations (28) 
and (29). The inner-most error bars represent the standard deviation of relative differences defined in 
equation (3 1). The standard deviation error bars indicate where most of the relative differences occur, while 
the minimum/maximum error bars indicate the extreme values of relative differences. The absolute value of 
all relative differences is clearly within 5%, and most relative differences are within 2%. This indicates that 
the simpler nearest-neighbor interpolation approach gives results that are very similar to the more 
complicated and computationally expensive log-log cubic Lagrange method. It should also be noted that the 
relative differences shown in Figures 2-4 translate into relative differences for dose equivalent that are 
bounded by 0.3%. It will be shown later, that the increased speed provided by the nearest neighbor 
approach justifies this negligible loss in accuracy. 



Figure 2. Relative difference between fluxes computed by the original and modified heavy ion transport 
algorithm as a function of kinetic energy in aluminum and water exposed to solar minimum GCR 
conditions. Light and dark error bars represent extreme and standard deviation of relative differences, 

respectively. 
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Figure 3. Relative difference between fluxes computed by the original and modified heavy ion transport 
algorithm as a function of shielding thickness in aluminum and water exposed to solar minimum GCR 
conditions. Light and dark error bars represent extreme and standard deviation of relative differences, 

respectively. 



Figure 4. Relative difference between fluxes computed by the original and modified heavy ion transport 
algorithm as a function of charge in aluminum and water exposed to solar minimum GCR conditions. Light 
and dark error bars represent extreme and standard deviation of relative differences, respectively. 


6. Convergence Test 

A computational algorithm is said to converge if the numerical solutions reach an asymptotic limit 
as the discretization parameters approach zero. HZETRN has two discretization parameters: the step size, /?, 
and the number of energy grid points, N. In order to show convergence and quantify discretization error, 
discretization parameters are refined several times and the differences between the various solutions are 
compared. This process is referred to as a convergence test. 

In this section, convergence tests are applied to the original and modified heavy ion transport 
algorithms to show convergence and quantify overall discretization error. Five different energy grid sizes 
are considered ranging from N= 100 to N = 502, along with eleven different step-sizes ranging from h = 
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0.5 g/cm" to h = 2" u g/cm 2 . These discretization parameters are used in both the original and modified 
heavy ion transport algorithms to transport a solar minimum GCR environment through 100 g/cm 2 of 
aluminum and water. The solar minimum GCR environment was calculated with the Badhwar-O'Neill 2010 
GCR model [O'Neill 2010] using a solar modulation parameter of 475 MV. The same light particle 
transport algorithm was coupled to both the original and new heavy ion transport algorithm. Dose 
equivalent is computed at a number of depths in the slab and is used as the quantity of interest for the error 
analysis. 

Dose equivalent as a function of depth in aluminum and water exposed to the solar minimum GCR 
environment is given in Figure 5 for both algorithms and all of the discretization parameters. The data in 
Figure 5 show that both algorithms with the various discretization parameters give very similar results over 
the thicknesses considered. The spread in dose equivalent values at 100 g/cm" in aluminum and water is 
1.2% and 1.4%, respectively. This variation includes both the differences associated with the original and 
new heavy ion transport algorithm as well as the inherent discretization errors. 



Figure 5. Dose equivalent versus depth in aluminum (left pane) and water (right pane) exposed to solar 
minimum GCR conditions computed with the original and modified heavy ion transport algorithms and all 

discretization parameters. 



Figure 6. Total discretization error for dose equivalent as a function of depth in aluminum (left pane) and 
water (right pane) exposed to a solar minimum GCR environment. 
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The total discretization error associated with the original and modified heavy ion transport 
algorithms for aluminum and water targets exposed to the solar minimum GCR environment is shown in 
Figure 6. For each algorithm, the discretization error is obtained by computing the relative percent 
difference between the converged results obtained with fine discretization parameters ( h = 2" 1 1 g/cm 2 and N 
= 502) and the solution obtained with typical discretization parameters (/? = 0.5 g/cm 2 and N = 100). A 
negative discretization error indicates that the results obtained with the typical discretization parameters 
underestimate the converged result. It is clear that the modified heavy ion transport algorithm results in a 
negligible increase in discretization error in both aluminum and water targets as thick as 100 g/cm 2 . It will 
be shown in the next section that this minor loss of accuracy is justified by the significant reduction in code 
run time. 

7. Run-time Comparisons 

The original and modified heavy ion transport algorithms are used to transport the same solar 
minimum GCR boundary condition through three different shielding configurations typically used with 
F1ZETRN. The solar minimum GCR environment was calculated with the Badhwar-O'Neill 2010 GCR 
model [O'Neill 2010] using a solar modulation parameter of 475 MV. These run-time tests were performed 
on a Linux computer with 8 quad-core AMD Opteron 8378 processors and 256 GB of RAM. Each test was 
run 25 times with the original and new algorithms, and the reported run-times are the average of the 25 
separate executions. Results are given in Table 1. The column labeled "1 layer" corresponds to a shielding 
configuration made of a single layer of aluminum. This type of calculation might be done to analyze the 
radiation protection properties of a slab of some material or to generate a simplified interpolation database 
for vehicle analysis. The column labeled "2 layer" corresponds to a shielding configuration made of an 
aluminum shield followed by a tissue target. Fluxes are stored on a two-dimensional grid of points covering 
every combination of first and second layer thicknesses. This simulation is typically used to model a 
shielding structure (aluminum) surrounding a human phantom (tissue). The column labeled "3 layer" 
corresponds to an aluminum shield followed by a polyethylene shield and a tissue target. This simulation is 
typically used to model a shielding structure (aluminum) with a secondary shield (polyethylene) 
surrounding a human phantom (tissue). In all three cases, the maximum thickness of any layer is 100 g/cm 2 , 
and each layer has 11 spatial grid points spaced nearly logarithmically at which fluxes are saved and 
printed to a file. Results are given in Table 1. 

Table 1. Run-time (se conds) comparisons between original and new heavy ion transpo rt algorithms. 



1 layer 

2 layer 

3 layer 

Original 

9.2 

96.1 

1042.3 

New 

2.9 

20.1 

212.4 

Factor 

3.2 

4.8 

4.9 


The overhead computational cost of reading in cross section databases and writing out flux files is the same 
in both codes and is included in the total run-time. Consequently, the ratio between the original and new 
run times is not constant. Two additional 1 -layer tests were run, but the slab thickness was increased to 
1000 g/cm 2 and 5000 g/cm 2 , and the fluxes were only printed out at the front and back boundaries of the 
slab to minimize the run-time dedicated to file writing. These tests help minimize the impact of reading and 
writing data on total run-time and more clearly show the increased speed associated with the new 
algorithm. It was found that the new code was 5.5x and 5.8x faster than the old code in the 1000 g/cm 2 and 
5000 g/cm 2 slabs, respectively. 

8. Conclusions 

In this work, a modification to the heavy ion transport algorithm in HZETRN was shown to be 
approximately 5x faster than the original algorithm, and convergence tests showed that the loss in 
numerical accuracy was negligible (< 1%). This increase in overall code run time for GCR calculations will 
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be beneficial for future analysis including risk assessment, ray-by-ray calculations, and geometries with 
large shielding thicknesses such as the International Space Station and martian atmosphere. 
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